home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_6 / a6_e.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  4.5 KB  |  179 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 6.e (Error Analysis for Numerical Differentiation).
  9. % Section    6.1,    Approximating the Derivative, Page 326
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program investigates the error E(f,h) when
  13. % computing numerical approximation for f'(x).  Two central
  14. % Two central difference formulas are explored:
  15.  
  16. %            f(x+h) - f(x-h)
  17. %  f'(x)  =  ---------------  +  E(f,h)
  18. %                  2h
  19.  
  20. %            - f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)
  21. %  f'(x)  =  ---------------------------------------  +  E(f,h)
  22. %                              12h
  23.  
  24. %    The function f(x) is stored in the M-file  f.m.
  25. % function z = f(x)
  26. % z = cos(x);
  27.  
  28. delete f.m
  29. diary f.m; disp('function z = f(x)');...
  30.            disp('z = cos(x);');...
  31. diary off;
  32.  
  33. f(0); % Remark. f.m is used for Algorithm 6.e
  34.  
  35. pause % Press any key to continue.
  36.  
  37. clc;
  38. %                f(x+h) - f(x-h)
  39. % For  f'(x)  ≈  ---------------
  40. %                      2h
  41.  
  42. %    The error bounds is known to have the general form:
  43.  
  44. %         |E(f,h)| ≤ eps/h + m h^2/6
  45.  
  46. % and the optimum step size is  
  47.  
  48. %         h = (3 eps/m)^(1/3)
  49.  
  50. % where  eps is machine epsilon (round off).
  51.  
  52. pause % Press any key to continue. 
  53.  
  54. clc;
  55. % Place the differentiation point in  x 
  56.  
  57. % Place a bound for   |f'''(c)|   in  m
  58.  
  59. x = 0.8;
  60.  
  61. m = 1;
  62.  
  63. hopt = (3*eps/m)^(1/3);
  64.  
  65. emin = eps/hopt + m*hopt^2/6;
  66.  
  67. pause    % Press any key to graph the error bound.
  68.  
  69. clc;
  70. a = hopt/10000;
  71. b = 3.44*hopt; %b = 0.00003;
  72. c = 0;
  73. d = 5.25*emin; %d = 0.0000000002;
  74. h  = (b-a)/150;
  75. H = a:h:b;
  76. B = eps./H + m*H.^2/6;
  77. E = abs(- sin(x) - (f(x+H) - f(x-H))./(2*H));
  78. axis([a b c d]);...
  79. plot(H,E,'-r',H,B,'-g');...
  80. hold on;...
  81. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  82. xlabel('h');...
  83. ylabel('e');...
  84. title('|E(f,h)| ≤ eps/h + m h^2/6');...
  85. grid;...
  86. axis;...
  87. hold off;...
  88. shg; pause    % Press any key to continue.
  89.  
  90. clc;
  91. Mx1 = 'Graph of the error bound  B = eps/h + m h^2/6,';
  92. Mx2 = 'and error in the approx.  E = |E(f,h)|.';
  93. Mx3 = 'The two curves verify the relationship:';
  94. Mx4 = '     |E(f,h)| ≤ eps/h + m h^2/6';
  95. Mx5 = 'The optimum step size from the theory is:';
  96. Mx6 = '     hopt = (3 eps/m)^(1/3) = ';
  97. Mx7 = [Mx6,num2str(hopt)];
  98. Mx8 = 'The minimum error from the theory is:';
  99. Mx9 = '     emin = eps/hopt + m hopt^2/6 = ';
  100. Mx10 = [Mx9,num2str(emin)];
  101. clc,echo off,diary output,...
  102. disp(''),disp(Mx1),disp(''),disp(Mx2),...
  103. disp(''),disp(Mx3),disp(''),disp(Mx4),...
  104. disp(''),disp(Mx5),disp(''),disp(Mx7),...
  105. disp(''),disp(Mx8),disp(''),disp(Mx10),diary off,echo on
  106.  
  107. pause % Press any key to continue.
  108.  
  109. clear; clc;
  110. %                - f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)
  111. % for  f'(x)  ≈  ---------------------------------------
  112. %                                  12h
  113.  
  114. %    The error bounds is known to have the general form:
  115.  
  116. %         |E(f,h)| ≤ 3/2 eps/h + m h^4/30   and the
  117.  
  118. % Optimum step size  
  119.  
  120. %         h = (45/4 eps/m)^(1/5)
  121.  
  122. % where  eps is machine epsilon (round off).
  123.  
  124. pause % Press any key to continue.
  125.  
  126. clc;
  127. % Place the differentiation point in  x 
  128.  
  129. % Place a bound for  |f'''''(c)|  in  m
  130.  
  131. x = 0.8;
  132.  
  133. m = 1;
  134.  
  135. hopt = ((45/4)*eps/m)^(1/5);
  136.  
  137. emin = (3/2)*eps/hopt + m*hopt^4/30;
  138.  
  139. pause    % Press any key to graph the error bound.
  140.  
  141. clc;
  142. a = hopt/10000;
  143. b = 2.51*hopt; %b = 0.003;
  144. c = 0;
  145. d = 5.78*emin; %d = 0.000000000002;
  146. h  = (b-a)/150;
  147. H = a:h:b;
  148. B = (3/2)*eps./H + m*H.^4/30;
  149. E = abs(- sin(x)-(-f(x+2*H)+8*f(x+H)- 8*f(x-H)+ f(x-2*H))./(12*H));
  150. axis([a b c d]);...
  151. plot(H,E,'-r',H,B,'-g');...
  152. hold on;...
  153. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  154. xlabel('h');...
  155. ylabel('e');...
  156. title('|E(f,h)| ≤ 3/2 eps/h + m h^4/30');...
  157. grid;...
  158. axis;...
  159. hold off;...
  160. shg; pause    % Press any key to continue.
  161.  
  162. clc;
  163. Mx1 = 'Graph of the error bound  B = 3/2 eps/h + m h^4/30,';
  164. Mx2 = 'and error in the approx.  E = |E(f,h)|.';
  165. Mx3 = 'The two curves verify the relationship:';
  166. Mx4 = '     |E(f,h)| ≤ 3/2 eps/h + m h^4/30';
  167. Mx5 = 'The optimum step size from the theory is:';
  168. Mx6 = '     hopt = (45/4 eps/m)^(1/5) = ';
  169. Mx7 = [Mx6,num2str(hopt)];
  170. Mx8 = 'The minimum error from the theory is:';
  171. Mx9 = '     emin = 3/2 eps/hopt + m hopt^4/30 = ';
  172. Mx10 = [Mx9,num2str(emin)];
  173. clc,echo off,diary output,...
  174. disp(''),disp(Mx1),disp(''),disp(Mx2),...
  175. disp(''),disp(Mx3),disp(''),disp(Mx4),...
  176. disp(''),disp(Mx5),disp(''),disp(Mx7),...
  177. disp(''),disp(Mx8),disp(''),disp(Mx10),diary off,echo on
  178.  
  179.